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ABSTRACT 

Determining the physical characteristics of a star is an inverse problem consisting in estimat- 
ing the parameters of models for the stellar structure and evolution, knowing certain observ- 
able quantities. We use a Bayesian approach to solve this problem for a Cen A, which allows 
us to incorporate prior information on the parameters to be estimated, in order to better con- 
strain the problem. Our strategy is based on the use of a Markov Chain Monte Carlo (MCMC) 
algorithm to estimate the posterior probability densities of the stellar parameters: mass, age, 
initial chemical composition,. . . We use the stellar evolutionary code ASTEC to model the star. 
To constrain this model both seismic and non-seismic observations were considered. Several 
different strategies were tested to fit these values, either using two or five free parameters in 
ASTEC. We are thus able to show evidence that MCMC methods become efficient with re- 
spect to more classical grid-based strategies when the number of parameters increases. The 
results of our MCMC algorithm allow us to derive estimates for the stellar parameters and 
robust uncertainties thanks to the statistical analysis of the posterior probability densities. We 
are also able to compute odds for the presence of a convective core in a Cen A. When using 
core-sensitive seismic observational constraints, these can raise above ~40%. The compari- 
son of results to previous studies also indicates that these seismic constraints are of critical 
importance for our knowledge of the structure of this star. 
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1 INTRODUCTION 

In stellar physics, a Cen A holds a particular place. Located at 
1.3 pc from the Sun, it is the closest star to the Solar System. Its 
apparent magnitude is V = -0.01, making it a privileged target 
to obtain precise and accurate data from spectroscopy, photome- 
try, astrometry or interferometry. Having a G2V spectral type, it is 
physically close to the Sun and is thus of importance for physicists 
interested in the processes and mechanisms governing the structure 
and evolution of solar analogs. Finally, it is part of a triple sys- 
tem and forms a close, visual, binary with a Cen B. This allows to 
directly measure the mass of the star, which is an information of 
critical importance in stellar physics. 

For all these reasons, a Cen A has been the focus of many ob- 
servational campaigns and theoretical studies. The former have al- 
lowed to obtain very precise data for the star. We possess measure- 
ments of its atmospheric parameters (effective temperature, surface 
gravity, chemical abundanc es; |Chmielewski et al.|1992{ Neuforge- 
Verheecke & Magain 1997), luminosity (Soderhjelm 1999), radius 
jKe^ella^et^lJ|2003J, oscillation frequencies ([Bouchy & Carrier| 



constrain theoretical models for the stellar structure and evolution, 
leading to a significant improvement of our knowledge of the phys- 
ical properties of a Cen A, its mass, age, initial chemical composi- 
tion. 

Despite all these advances, it is obvious that there is no per- 
fectly unified picture yet of the internal structure of a Cen A, nei- 
ther of its evolutionary state. Even though its mass is well-known, 
tiiere exist discrepancies in the ages found in the literature |Tur-| 



2002, Bedding et al. 2004),... All these have been combined to « 
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cotte & Christensen-Dalsgaard 1998 ; Guenther & Demarque 2000; 
Morel et al.|2000 ; Eggenbe rger et al.|2004| |Miglio & Mon talban 
2005). It indeed depends very much on the considered observa- 
tional constraints and the particular physics included in the models 
that are used to reproduce these quantities. 

Furthermore, as it turned out, a Cen A is in a very peculiar 
state with regard to stellar structure and evolution. Indeed, stars 
with mass >1.1 M Q can be convectively unstable in their cores. 
This is very important because it affects the energy transport in the 
star and thus its entire evolution. However, it is a difficult task to 
identify such small convective cores. First, they might not influ- 
ence much the stellar observables (with the notable exception of its 
oscillation frequencies) and could remain undetected. Second, the 
physics of convection governing these structures is poorly known, 
and it is likely that the models used to describe a Cen A are in- 
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correct with respect to the treatment of its inner convection. This is 
reflected by the fact that various studies have led to different claims 



Regarding the existence of a convective core in a Cen A 


Thevenin 


et al.|2002||Thoul et al.|2003 pggenberger et al.|2004| 


Miglio & 


Montalban 2005 , de Meulenaer et al. 2010). 



A further problem appears when one needs to quantify pre- 
cisely the uncertainties on the physical stellar parameters estimated 
relatively to the uncertainties on the observations. This is not a 
problem specific to a Cen A: it has always been challenging in 
stellar physics to obtain robust and credible uncertainties on the pa- 
rameters estimated from modelling. This is mostly due to the sub- 
stantial computational cost of the stellar evolutionary codes, which 
is often incompatible with the intensive use of the models required 
by statistical methodologies. It does not come as a surprise then 
that most of the efforts in stellar physics have long been directed 
towards improving the codes themselves, either numerically or by 
including more accurate physics, rather than developing tools for 
parameter estimation. 

This situation has progressively changed during the past 
decades. Many studies have been published providing extensive 
grids of models, spanning large ranges of values of physical pa- 
rameters such as the stellar age, mass or initial chemical compo- 
sition jScharler et al.|1992| |Girardi et a"L]|2000[ |Pietrinferni et al.| 
|2004||Quirion, Christensen-D alsgaard & Arentoft 2010). These can 
help to proceed to statistical studies without having to compute too 
many models at once. Simultaneously, there has been a growing 
interest for methodologies for parameter estimation, which have 
been applied to stellar physics. Among them, optimization meth- 
ods have been introduced first, in some cases great care have been 
used to combine this procedure with uncertainty analysis such as 
SVD decomposition (Brown et al. 1994). The methods adopted 
for optimization are many, from straightforward grid analysis to 
solving the non-linear least-square problem using a Levenberg- 
Marquardt algorithm (Miglio & Montalban 2005) or stochastic 
sampling via genetic algorithms ( |Met calfe 2003 ). Other studies 
have focused on the application of Bayesian statistics with some 
success. Here again, the methodologies differ, some using grid- 
based strategies (J0rgensen & Lindegren 2005 , d a Silva et al.|2006| 
|Takeda et al.|2007) other Markov Chain Monte Carlo (MCMC) al- 
gorithms (Bazot, Bourguignon & Christensen-Dalsgaard 2008). 

In this article, we present a Bayesian modelling of a Cen A 
based on such an MCMC methodology. In Sect. [2] defining our 
inverse problem, we briefly recall the essentials of the Bayesian ap- 
proach. We also give a brief account on the MCMC theory and 
present the algorithm we used to solve this inverse problem. In 
Sect. [3] and [4] we present the physical context to which these tools 
are to be applied. We review briefly the characteristics of our stel- 
lar evolution code and give an overview of the current observational 
knowledge on a Cen A. In order to assess the potential gains one 
can expect from the use of MCMC algorithms with respect to grid- 
based strategies we give two solutions to our inverse problem. In the 
first one we let two stellar parameters vary and make assumptions 
to fix three others. In the second one, we estimate all five. We re- 
port the results from both numerical experiments in Sect. [5] In par- 
ticular, we concentrate on two points. One is physical: what is the 
general picture of the structural and evolutionary state of a Cen A? 
The other is numerical: how efficiently were these results obtained? 
From the methodological point of view, our results show the inter- 
est of MCMC approaches for such estimation problems, when no 
explicit mapping exists between observables and the parameters to 
be estimated. In particular, compared to grid-based computations, 
MCMC reveals better efficiency by concentrating most of the com- 



putational effort on the regions of interest. They also naturally pro- 
vide crucial information about uncertainties associated to the esti- 
mated values. We briefly compare our results to previous studies in 
Sect. [531 



2 THE BAYESIAN APPROACH 
2.1 A difficult inverse problem 

Stellar models include several parameters that have to be estimated 
by matching observational data, which represents a typical infer- 
ence problem. Stellar observable quantities usually considered are 
the luminosity L, the effective temperature r cff , the surface metal- 
licity Z, the radius R and the frequencies of stellar pulsation modes 
in a given range of radial orders n and angular degrees /, say (v„,/). 
a Cen A is one of the very important cases in stellar physics for 
which all of these measurements are available to high precision. 
From such observations, one usually wants to estimate stellar pa- 
rameters such as the mass of the star, M, its age, t t , the mixing- 
length parameter, a, and the initial mass-fractions of helium and 
heavy elements, respectively Y and ZoQThis is by no mean an ex- 
haustive set of parameters, in particular if non-standard physics are 
to be included. 

Unfortunately, there is no closed-form equation linking stel- 
lar parameters together with observables, which could be inverted. 
For a given set of interest stellar parameters, corresponding val- 
ues of L, r cff , Z, R and (v n/ ) can be computed by stellar evo- 
lution codes (which may also require additional parameters). Let 
6 = (9\, . . . ,9 K ) collect the unknown stellar parameters and X = 
(Xi, . . . ,X N ) collect the observed quantities. We shall call the pa- 
rameter space, P, and the observable space, O, the spaces contain- 
ing respectively and X. We have dim(P) = K. 

We shall now consider a mapping from P to O and write 

X = y(0\ 

where elements in X are the values predicted by the stellar code for 
these quantities, denoted here by 5?. We consider that the observa- 
tions are related to X by 

X = X + e, 

where e is an unknown error term accounting for instrumental er- 
rors on the measurements and possible model errors between X 
and y{6). The problem addressed here concerns the estimation of 
6 from 

X = 3?(0) + e, (1) 

which defines an inverse problem (Idier 2008). A classical ap- 
proach to solve this problem consists in minimizing an appropri- 
ate distance, in O, between X and X, such as the usual x 2 error 
||X-^(0)|| 2 = 2" =1 (X„-X„) 2 .Thisisahardtaskinourcase.First, 
SP{&) is not an explicit function, so no analytical formulation can 
be exploited by efficient optimization methods (for example, no ex- 
pression of the derivatives of with respect to each variable in 
6 can be used for optimization). Second, stellar models computed 
by make use of highly non-linear functions, so that classical 
optimization methods may be stacked in local minima or yield de- 
generate solutions. The complexity of the mapping from the space 
of parameters to the space of observables has already been pointed 

1 Or Xq, the initial hydrogen fraction, remembering that the relation Xq + 
Y + Zq = 1 holds. 
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5 OUt ( Brow n et al.p994) [j0r gensen & Lindegren|2005| > and repre- 
sents one of the major challenges in stellar-parameter estimation. It 
has to be noted that obtaining S?(ff) for a given set of parameters 
is generally computationally expensive. Therefore, computing the 
model on a grid of discrete values of the parameters may be ex- 

10 tremely time-consuming if one wants to explore the whole parame- 
ter space with high precision, especially as its dimension increases. 

We consider the use of Markov Chain Monte Carlo meth- 
ods (MCMC) as an alternative to classical optimization methods, 
where model computations are concentrated around the interest re- 

15 gions of the parameter space, and are in that sense "optimized". 
Basically, the distance between observations and their prediction 
by is linked to a probability distribution, and MCMC algo- 

rithms aim at approximating it by generating random samples dis- 
tributed according to such distribution. Then, they allow one to ob- 

20 tain not only an estimate of the interest parameters, but also to asso- 
ciate characterizations such as uncertainties or higher-dimensional 
statistics. MCMC are also naturally oriented toward a Bayesian set- 
ting, in which prior information on the unknown parameters is also 
taken into account by means of probability distributions. Further- 

25 more, being stochastic in nature, they might help to avoid confine- 
ment in local minima. 



2.2 Probabilistic framework 

2.2.1 Likelihood function 

Equation {T} has a statistical interpretation, where the perturbation 
30 term e can be seen as the realization of a random process. We sup- 
pose hereafter that the error e n on each observation X„ is normally 
distributed, with variance cr 2 . and that the errors affecting the differ- 
ent observations are independent one from each other. In our case, 
we assume that cr„ corresponds to the uncertainty on X„, which is 
35 determined by the observation. 

Then, the probability density of all perturbation terms taken 
together is 



(2tt) 



- exp 



which defines the likelihood function of data X given a model pa- 
rameterized by 0, using l[TJ 
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Then, maximizing the likelihood l|2j is equivalent to minimizing 
the usual least-squares term 



mm^ 2 (6l) = 



(3) 



2008 ), by considering such information through a prior probability 
45 distribution on 0, say n(0). 

Suppose, for example, that the age is bounded between t mm 



2.2.2 Incorporating prior information 

Suppose one has some - even weak - prior information about the 
40 unknown parameters 0, such as an acceptable range of possible val- 
ues or previous measurements with associated uncertainties. One 
may want to take it into account in order to constrain the solution 
of {3}. Bayesian statistics provide a natural way to do this ( |Idier| 
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,(f*), 



(4) 



where is the indicator function of the domain S>, which 

equals 1 if u 6 & and otherwise. That is, n(t t ) is a uniform distri- 
bution on Or suppose that the mass M has been previously mea- 
sured at M cst , with uncertainty cr Mca . Then, one may use a Gaussian 
prior distribution n(M), centred at M est and with standard deviation 
o~M ca (or 2cr Mca , 2>o- Mcsx for a weaker prior distribution). 

Once such priors are defined, can be estimated from the pos- 
terior probability distribution (PPD), which reads, by Bayes' theo- 
rem 

y(fl;XW9) 



n(0\X) : 



m(X) 



(5) 



n(0\X) represents the probability distribution of all unknown pa- 
rameters constrained by observed values of X and by prior distri- 
butions. In the latter equation, m(K) is a normalization constant en- 
55 suring that the integral of n(0\X) over the entire space of parameters 
is 1. 



2.2.3 Exploiting the posterior distribution 

Suppose we are able to access the whole posterior probability dis- 
tribution n(0\X). It contains all statistical information available, 
given a set of related observed quantities and (possibly) under the 
given prior assumptions. If no prior distribution is used, then n(0\X) 
identifies with the likelihood function Jzf(0;X). We are then theo- 
retically able to exploit such a distribution in order to derive esti- 
mates of interest parameters 0. In this paper, we focus on the two 
most common ones in Bayesian analysis. The most intuitive one is 
the Maximum A Posteriori (MAP) estimate, which is the set of pa- 
rameters that maximize the PPD (that is, "the most probable value" 
of given observations X) 

0map = ar g max n(0\X). (6) 

The Posterior Mean (PM) estimate corresponds to the expected 
value of under distribution {5} 



0p 



E[0|X] 



0n(0\X)dO. 



(7) 



It minimizes the mean square error E [(#(X) - 0) 2 1. 

The latter one is also often used in Bayesian estimation since, 
knowing n{0\X), one can associate variances, that is, uncertainties, 
to the estimated values. 

In practice, of course, we do not have access to the whole 
posterior probability distribution and computation of the former 
estimates is a hard task. Optimization of n(0\X) is at least as dif- 
ficult as optimization of likelihood |2j, for the reasons explained 
in Sect. |2.1| tt(0\X) does not have an exploitable closed-form ex- 
pression. Computing 0pm requires an integration over all possible 
values of and is also impossible analytically. 

2.3 Markov Chain Monte-Carlo algorithm 



MCMC methods (for a general review see |Robert & Casella|2005| > 
are powerful techniques allowing to approximate complex proba- 
bility densities. Roughly speaking, MCMC aims at generating sam- 
ples randomly distributed according to the interest distribution j3J, 
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and then approximating such distributions by histograms. MCMC 
methods are powerful in the sense that they are able to generate 
random samples distributed (at least, asymptotically) according to 
a target distribution only by evaluating this distribution at a certain 
number of parameter values. 

In our case, suppose we have T samples . . . , (f'\ . . . , (T> 
distributed according to n(0\X) - we use the notation (,} ~ n(6\X). 
Then, the PM estimate l|7j of all parameters in can be approxi- 
mated by 



Opm 



and the corresponding unbiased variance estimator by 



1 T - 
o~ 2 (#*,pm) - j, _ r ^ - #/S,pm) 



(8) 



(9) 



with 6k the k-th component of 0. We can also approximate the 
MAP estimate by taking, among all the set of parameters for 
which 7r(# (,) |X) is the highest (note that the corresponding values 
are naturally obtained as a by-product of the MCMC algorithm). 

MCMC methods work by generating a Markov Chain, that is, 
a sequence . . . ,# (r) such that the rule for generating de- 
pends only on the previous value O'-'^'K The Monte-Carlo refer- 
ence obviously concerns the random nature of the process. Giving 
a full account of the theory of MCMC is obviously far beyond the 
scope of this paper. In the following, we will only comment on 
some essential practical aspects necessary for the comprehension 
of the subsequent parameter estimation. This approach was first 
suggested in the framework of ste llar modelling by Bazo t, Bour-| 
guignon & Christensen-Dalsgaard ( 2008 1 



Algorithm [A.l] Metropolis-Hastings algorithm 

For t = 1, choose an initial value Then, at iteration 
V. 

1. Generate: 0' ~ 4(0|0 (,) ) 

2. Select: 

„_„ \o" with probability p(0 (,) , 0*), 
G — 



where 



(T with probability 1 - p{&'\ 0* ), 



p(0 ',0 ) = mmi - : — — , 1 



;r(6HX) q(fiW>) 



2.3.2 Instrumental distribution 

A key point conditioning the efficiency of the MH algorithm is the 
choice of the instrumental distribution. First, it should be compu- 
tationally easy to draw random samples distributed according to q. 
Second, q should concentrate the proposed values 0" around the 
interest regions of tt(0|X). Finally, q must also explore the whole 
support of possible values in order not to miss plausible parame- 
ter sets. It should be noted that the first condition represents one 
of the major computational interests of using MCMC algorithms: 
it becomes possible to simulate complex distributions using only 
simpler ones. 

Trying to ac count for these issues and following former sim- 
ilar choices (e.g., Andrieu & Doucet|l999 we choose g(#*|fl (f) ) as 
a mixture of two Gaussian distributions centred at <f> and the uni- 
form distribution on the definition domain of all parameters. That 
is, for Step 1 in the MH algorithm, we generate, for all parame- 
ters with k = l,...,K 



2.3.1 The Metropolis-Hastings algorithm 

In this paper, we construct a Metropolis-Hastings (MH) algorithm 
(Metropolis 1953, Hastings 19701, which is the most generic class 
of MCMC methods. The generality of this method allows to esti- 
mate probability densities of all forms, with the only assumption 
that they can be evaluated at every point in the parameter space. 
The PPD expressed in {5} falls within this category. 

The Markov Chain is constructed through a two-step iterative 
procedure. Suppose that . . . , # (T) have been obtained. First, a 
trial value of the set of parameters, 0*, is generated according to 
an instrumental distribution q(0*\0 f ' } ), that may depend on the pre- 
vious value of the chain Then a selection process operates, 
where # <,+I) is set to 0* with some probability p, and to other- 
wise. The convergence of the distribution of {0 1 '*} toward the target 
distribution is then ensured by the analytical expression of p, see 
the algorithm description below. We give more details about moni- 
toring of the convergence of MCMC in Appendix [A] 

The choice of 0* at step 1 is detailed in § |2.3.2| In practice, 
accepting = 0* with probability p is performed by generat- 
ing u according to a uniform distribution in [0, 1], and # (,+ ' = 0* is 
accepted if u < p. 

Note that if q is a symmetrical distribution (that is, q(0 {,) \0") = 
q(6*\w'')), which will be the case here, then the ratio in probabil- 
ity p is nothing but the PPD ratio evaluated at 0* and Hence, 
parameters increasing the PPD are always accepted while parame- 
ters decreasing the PPD are randomly accepted, the probability of 
acceptance decreasing with the PPD ratio. 



k >H,2 



) with probability p\ , 
with probability p2, 
with probability pi . 



(10) 



where jV^ e is the Gaussian distribution with mean 

and variance , , truncated to the definition domain % k of 6t, and 
is the uniform distribution on S) St - see Sects 



5.1 



and 



5.2 



which we specify the p, and the acceptable definition domain 
for 0. 

The use of Gaussian distributions centred on the current value 
# (,) defines random walks in the space of parameters. Having two 
of them helps to perform local scans of the parameter space at dif- 
ferent scales. We added the uniform distribution to avoid possible 
trapping in local extrema: if the chain is being stuck in such an area, 
there is chance for the algorithm to propose, from time to time, val- 
ues of well outside this region. 



3 STELLAR MODELS 
3.1 Physics 

In Eq.[T]appears the quantity S^{0), describing the predicted values 
for the observable quantities. These values were computed using 
the Aarhus STellar Evolution Code (ASTEC). A thorough descrip- 
tion is given by Christensen-Dalsgaard ( 1 982 , 2008b ), we thus only 
briefly recall the main outline of the physics used in our models. 
The oscillation frequencies were computed with the adipls pack- 
age! Christensen-Dalsgaard 2008a). We restrict ourselves to the so- 
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called standard physics. This implies a spherical, non-rotating star. 
Microscopic diffusion, magnetic field and mixing processes other 
than convection (overshooting, rotational mixing,...) are neglected. 
10 J^his latter is described using the mixing-length theory ( |Bohm-| 
[Vitense|1958) that requires the setting of a mixing-length param- 
eter, a, defined as a = f„JH p , where £ m is the mixing-length and 
H p is the pressure scale-height.. 

The physics used in the stellar code is relatively basic. It is not 
15 the goal of this study to adopt the state-of-the-art physics for the 
stellar evolution code. The emphasis was rather put on the need to 
use a fast, robust and already tested version of the code. This is the 
jeason why we used the EFF tab les for the equation of state (Eggle- 
ton, Faulkner & Flannery 1973) instead of the latest OPAL tables 
20 (Rogers, Swenson & Iglesias 1996 ; Rogers & Nay fonov 2002), not- 
ing however that they agree "reasonably well". It also appears from 
Miglio & Montalban ( 2005 ) that, in the framework of the modelling 
of a Cen A, using EFF does not give results departing significantly 
from those using OPAL. Other input physics included the NACRE 
25 nucle ar reaction rates jAngulo et al."fT9 99 ) and OPAL opacity ta- 
bles Pglesias &"R ogers 1995}. The relative abundances of heavy 
elements are those prescribed by |Grevesse, Noels & SauvaT|jl993| 
GNS abundances). 

This latter point should be briefly commented since these 
30 abundances are, at the moment, subject to debate. Using 3- 
dimensional simulations, Asplund et al. (2004) have derived new 
solar abundances (AGS abundances). The most remarkable fea- 
tures are the low values found for C, N, O. Abundances of other 
heavy elements such as Na, Mg or Fe were also revised to lower 
35 values. The overall impact is a reduction of ~ 30% of the estimated 
solar metallicity. However, including these values in solar models 
has led to sound-speed profiles in disagreement with helioseismol- 
pgy (e.g. | Turck-Chieze et al.|2004"[|Montalban et al.|20"04"[|Antia &] 



Table 1. Observational constraints on a Cen A. The last column specifies 
the reference for the selected quantity. 



[PS. . 

|Basu||2"005) . Although several attempts have been made, some of 
40 them including non-standards physics (e.g., |Guzik, Watson & Cox] 
2005 ; Castro, Vauclair & Richard 2007 ), no satisfactory agreement 
between theoretical models and seismology could be reached. On 
the contrary, models using the GNS abundances, computed using 
1-dimensional models, can reproduce with a remarkable precision 
45 the solar sound-speed profiles. It should be noted that recent studies 
(e.g., |Caffau et al. 2008 , As plund et al.|2 009 ) have led to upward re- 
visions of these abundance. Discriminating between AGS and GNS 
abundances is beyond the scope of this study. We thus keep on us- 
ing the GNS heavy-element abundances in our models of a Cen A. 
so This rather conservative course of action is motivated by the fact 
that we focus on testing a new technique for statistical inference 
of stellar parameters. That should be done with input physics that 
has already led to coherent results. At the very least, the inferred 
parameters for a Cen A could be compared to other stellar models 
55 that also used GNS abundances. 

Given these physical prescriptions, it is important to note that 
the "natural" parameter space is such that dim(P) = 5, with the 
typical associated vector 6 = (M, t+,a, Zq, Xq). 



4 OBSERVATIONAL CONSTRAINTS 

60 In the following we review briefly the values we selected for the 
different observables for a Cen A, their uncertainties and their po- 
tential bias. We present both the non-seismic and the seismic ob- 
servations. Table [T] sums up the former ones. For the oscillation 
frequencies, the reader is referred to Bazot lit"aL"] ( |2007| >. 



Observable 


Value 




Reference 


T e g 


5810 ± 50 K 


H 


;genberger et al. 


{2004 


L/L 


1.522 ±0.030 




jgenberger et al. 


(2004 


R/Re 


1.224 ±0.003 


Kervella et al. 


(2003 


Z/X 


0.039 ± 0.006 




|Thoul et al. 


(2003 


<Av> 


105.9 ±0.1 fiUz 




Bazot et al. 


{2007 


(Sv) 


6.9 ± 0.4 pHz 




|Bazot et al. 


(2007 


M/M e 


1.105 ±0.007 




|Pourbaix et al. 


(2002 



65 4.1 Spectro-photometric observables 

Several spectroscopic determinations of the effective temperature 
of a Cen A have been carried out, whether by fitting the wings of 
the H„ line or the Fel lines (Neuforge-Verheecke & Magain 1997 ). 
However some discrepancies appear between the value used in the 
70 modelling by |Morel etal~| ( |2000) , 5790 ± 30 K obtained by fitting 



the wing of the H„ line, and the value given by Neuforge-Verheecke 
jfe Magain| ([1997), 5830 ± 30 K. To solve this problem, [Eggen- 
berger et al. (2004) have 



2004) have adopted an average value of these two 



measurements and lcr error bars encompassing their uncertainties; 
the same value was used by Miglio & Montalban (2005). Since 
it can help to compare our results with these works, we also use 
T eff = 5810 ± 50 K, even though reduced error bars of ±30 K could 
have helped reducing the uncertainties on the estimated parameters. 

Although the proximity of a Cen A should allow a precise 
measurement of its parallax, discrepancies also subsist between dif- 
ferent studies, affecting the subsequent values of the luminosity. 
We select the parallax determined by Soderhjelm ( 199 9). The cor- 
responding luminosity is L/L Q = 1.522 ± 0.030 ( Eggenberg er et ah] 
|2004[ ). Note that this value of the parallax has also been used to 
derive the radius and the mass we describe below. 



4.2 Radius 



The radius of a Cen A has been measured by Kervella et al. ( 2003 ) 
with VINCI, the VLTI Commissioning Instrument, which operates 
in the K-band (2.0-2.4 /im). They measured an angular diameter, 
after correcting for limb-darkening effects, d = 8.5 1 1 ± 0.020 mas. 
Combined with the parallax derived by Soderhjelm ( 1999), it led 
to R = 1.224 ± 0.003 R Q . For the reference solar radius value, 
jj^ye used the value mea sured from helioseismology by Br own &| 
[Christensen-Dalsgaard| \199S\ , R e = 695.508 Mm. It is notewor- 
thy that the authors argued that this value, measured in the near- 
IR, can be compared to models computed with the stellar evolution 
code CESAM (Morel 1997) since the radius is there defined as the 
point where T = T s g , which is close to the temperature of the layer 
at which the continuum at 2.2 /im is formed. Since ASTEC uses 
the same atmospheric boundary conditions I Christensen-Dal sgaard| 
1982), we consider that this assumption still holds in our case. 

This radius value has been used in Eggenberger et al. (2004), 
along with seismic data from Bouchy & Carrier ( 2002 ), to constrain 
their model. They concluded that they cannot reproduce at the same 
time the seismic measurements and the radius within their lcr error 
bars. However, they noted that this could be done when consider- 
ing a 2cr interval. This confirmed the 1 ,3cr agreement Kervella et al. 
( |2003) found between the observed value of the radius and a pre- 
vious theor etical estimate using CORALIE frequencies ( |Thevenin"| 
[et al.|2002[ . |Miglio & Montalban| ( |2005| > did not use the observed 
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radius as a constraint for their models but always obtained values 
higher than the measured one. |Montalban & Miglio (2006), using 
this constraint in combination with seismic data from Bedding et al. 
(2004), were able to reproduce all observations within lcr. 



4.3 Seismic observations 



expanded the eigenfrequency equation to higher orders 



4.3.1 Individual frequencies 



The frequencies used in this paper are those from |Bazot et al.| 
( |2007} . As far as their accuracy is concerned, they can be con- 
sidered as reliable. The high SN ratio that can be achieved in the 
radial-velocity spectrum of a Cen A allows to find easily the most 
prominent modes, whose amplitude can reach up to 14 times the 
averagenoise level. A comparison between the studies byjBouchy 
& Carrier| ( [2002) ; |Bedding et al.| ( |2004) and |Bazot et aT1j2007|> of- 



fers a convincing empirical evidence that the modes are indeed real. 
The identifications presented in these studies are consistent (which 
of course does not prove them correct). 

Of much greater difficulty is the task of estimating confidence 
intervals on the estimated frequencies. The assumption made in 
|Bazot et al.| < [2007> may be viewed as the crudest one: considering 
that the modes are unresolved, the frequency resolution is chosen 
as the uncertainty on the frequency. The generic uncertainty given, 
1.3 /jHz, may be viewed as conservative. At the time our MCMC 
simulations were first performed, this was as good an estimate as 
one could find. We note however that methods are currently un- 
der development to allow better and more trustworthy estimation 
of the confidence intervals on the oscillation frequencies of solar- 
like stars (Brewer et al.|2007[|Stahn & Gizon|2008||Benomar et al.| 
|2010| [Handberg & Campante 2011) A recent discussion on this 
methods has shown that they are efficient even for ground-based 
gaped and irregularly sampled data ( |Bazot et al.|2012| >. While keep- 
ing in mind that room for potential improvement exists, we never- 
theless used the original uncertainties in the present study, since 
their values may at the very least represent a realistic upper bound. 
We should also note that a new seismic study of a Cen A using 
the previous data of |Bouchy & C arrier ( 2002 1 and Beddin g et al.| 
(20041 has been published by |de Meulena er et al.| ( |2010| > during 
our study. They have not been included here, but we note that their 
result partially agree with those of Bazot et al. (20071. 

Finally, we should note that fitting the individual frequencies 
is difficult because they are affected by surface effects that we 
cannot model. There have been tentatives to bypass this problem 
(Kjeldsen, Bedding & Christensen-Dalsgaard 2008) and thus allow 
to use full sets of frequencies as an observational constraints. These 
methods have gained in popularity and are often used in the frame- 
work of Kepler for the modelling of stars. However, these strategies 
are empirical and we chose not to use them in the present study. We 
instead rely on the frequency combinations that we present below. 
A possible development to be tested along these lines would con- 
sist in including in our estimation the Bayesian formalism given 
by Gruberbauer, Guenther & Kallinger (20121 that deals with this 
problem. This is left for future studies. 



4.3.2 Large separations 

The first asymptotic relation was given in the limit of low-degree, 
high-order p modes by Vandakur ov| ( |1967} . |Tassoul|(T9 80 ) further 
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and with R, the upper turning point of the mode and B a surface 
phase term. 

If we assume, to the lowest order of approximation, that the 
phase-related term quantity e s does not depend on the frequency, 
we immediately see that modes of consecutive order and same de- 
gree are equally spaced. The difference Av„j = v„j - v„_i_/ is called 
the large separation and may be considered roughly as a constant. 
This has been a useful tool in helio-and asteroseismology to iden- 
tify pulsation modes (e.g., Grec, Fossat & Pomerantz 1983 , Bouchy 
[& Carrier|2002} . 

The large separation is, according to Eq. l )12[ l, the inverse of 
twice the acoustic radius (i.e. the travelling time of an acoustic 
wave from the centre to the surface). We thus obtain a new global 
quantity that can be used to constrain stellar models. As noted by 
|Deubner"& Gough ( 1984), this integral depends mostly on the ex- 
ternal regions of the star. This can be easily seen if one considers 
that the sound speed in the stellar interior is a roughly monotonic, 
smoothly decreasing function of the radius, and that the ratio of its 
values at the centre and at the surface for a solar-like star is ~ 70. 
Therefore any acoustic wave will spend a much longer time closer 
to the surface than to the centre and will be predominantly sensitive 
to this region. 

It has been been noted that the large separation depends on 
the frequency. Allowing for a the surface-phase term to vary with 
frequency, in an equation similar to |TTJ, Vorontsov ( 1991 1 pointed 
out that it will control the overall behaviour of Avo when seen as 
a function of v. Once more, some problems arise when it comes 
to stellar models: as noted by Vorontsov & Zharkov ( 1989}, they 
usually provide poor values for this surface term and its frequency 
derivative, which is directly related to Av . 

From the point of view of stellar parameter estimation, this is 
nevertheless an interesting quantity. Besides the new (and relatively 
easy to extract) global constraint it provides, it is possible to com- 
pare the stellar parameters inferred using the large separations and 
other seismic indicators, obtained using different frequency combi- 
nations sensitive to deeper layers. This may in turn give us some 
insight on the modelling errors occurring below the surface. The 
quantity we used is the average large separation that we define as 
(Av) = Avo and compute as explained in |Bazot et al. (2007} and 
note that the estimate they give should read 105.9 ±0.1 /jHz and 
not 105.9 ± 0.3 yuHz; we therefore adopt the former value. 

Numerically, we computed the average large separation by 
performing a linear regression, for each degree / = 0, 1,2, of the 
relation v = v(n), only selecting in the code output the modes with 
orders corresponding to the values observed for a Cen A. This gives 
us three slopes, which are the average large separation for each de- 
gree. We then averaged over / to obtain the average large separation. 
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io 4.3.3 Small separations 

The second quantity used as a seismic indicator comes also 
straightforwardly from Eq. (TTV It immediately appears that modes 
separated by one one order and two degrees, v n j and v n ,-\,i+i, are de- 
generate at the first order. The difference between these frequencies 

15 is a second order quantity in 1 /v. It is called the small separation. 

Contrary to the large separations, the small ones depend 
mostly on the deep interior of the star. This can be seen from 
Eq. < | 1 3 1 > . The integral term behaves as c/r which is large at the 
centre and negligible near the surface. In-depth analysis by[Gough 

20 ( [T986) , [GabrieTl ( fT989l > and |Van Hoolst & Smeyers| ( p9lj > described 
much more rigorously this dependence. 

Small separations are obviously a very important tool for stel- 
lar parameter estimation using seismology. It is noteworthy that 
[Christensen-Dalsgaard] ( |1993[ > has shown that, in principle, the 

25 knowledge of the average values of the large and small separations 
is sufficient to characterize fully the evolutionary state of the star. 
We performed test estimation using the individual and the aver- 
age large separations. The main difference is that the small uncer- 
tainty on the latter provides a better constraint to our stellar model, 

30 therefore, we adopted the value provided by B azot et al"] |2007 ), 
(Sv) = 6.9 ± 0.4 /jHz. The average small separation in our models 
was simply computed by averaging our theoretical individual sep- 
arations over the mode order (selecting only the observed ones). 
- Finally, we als o tested the use of the ratio describe by |Rox-| 

35 [burgh & Vorontsov (1994), which should be completely free of sur- 
face effect^ 2 ] However, the large uncertainty on these quantities, in 
the range 22.4%-30.2%, made them poor constraints in the case of 
this a Cen A data set. 



55 5 RESULTS ON a Cen A 

5.1 MCMC simulations with dim(P) = 2 

5.1.1 Implementation details 

Our first set of MCMC simulations aimed at sampling the PPD in a 
parameter space of small dimension. This has some advantages, in 

60 particular for the comparison with procedures involving direct in- 
tegration of the normalization factor in equation Eq. <|5j from fixed 
grids of models (Sect |5.l!4| . It is indeed possible, in that case, to 
compute very dense grids of stellar models. We thus limit ourselves 
to a parameter space such as dim(P) = 2, with the stellar mass and 

65 age the only two free quantities. 

In order to fix the other parameters, we had to make a 
certain number of assumptions. The initial hydrogen fraction was 
considered roughly solar X = X 0o ~ 0.7. Combining with the 
[Fe/H] ratio for a Cen A, and because we neglected diffusion, 

70 this implies Z = 0.027. This of course is an oversimplification 
since we intentionally neglect the effect of the uncertainty on 
[Fe/H]. However, such a configuration may happen if one is, 
for instance, using grids of stellar models only sparsly sampling 
the initial metallicity parameter. The mixing-length parameter is 

75 somewhat more problematic since there is no fundamental physical 
justification for introducing this parameter. It can be related to the 
gradient of entropy in the considered convective zone, but there 
is no physical prescription that can be unequivocally applied to 
this quantity. It is possible to use a canonical value, calibrated 

so from the solar case. This is the approach we used here, setting 
or = a ~ 1.8. Some numerical simulations suggest that the value 
of the mixing-length parameter in the solar neighbourhood in the 
(M - /*) plane is close to constant (Trampedach 2007) . 



4.4 Using a mass measurement as a prior 

a Cen A being part of a close binary, its mass can be measured us- 
ing astrometry. Although its period is somewhat long, observations 
have been carried out throughout the last century, culminating with 
the precise measurements of Murdoch & Hearnshaw ( 1993 1, Endl 
[eTallpOOTT l and |Pourbaix et al.H2002) T using high-precision spec- 



rometers. This led to the derivation of a precise mass by Pourbaix 
iTaL] ( [2002l >, M obs = 1.105 ± 0.007 M (the error being give: 
lcr level). This mass was also used by Eggenberger et al. 



Migl io & Montalban|p005) and |Montalban & Miglio|p0Qi 



j2004] >, 

01. 



We used this value to obtain prior information on the mass. 
The assumption is made that the uncertainties are Gaussian and 
reflect truly the lcr deviation from the mean value, therefore 



7t(M) oc exp - 



(M - M obs y 



(14) 100 



It should be noted that Pourbaix et al. ( 2002 1 issued a warning con- 
so cerning the accuracy of this value. Even though the mass ratio be- 
tween a Cen A and a Cen B can be accurately constrained, the 
same is not true for their individual masses, because of the com- 
bined effects of the gravitational redshift and convective blueshift 
(this latter being especially difficult to correct for). 



2 Whereas there is still a contribution of the surface layers to the small 
separations. In practice, our test runs with this so-called Roxburgh's ratio 
on a Cen A gave us the same estimates as those using the small separations, 
albeit with much larger error bars. 



The set-up was similar for all our MCMC simulations. 
For each combination, X, of the observables we considered, we 
launched m = 3 independent Markov chains of length T = 50000. 
The use of parallel chains is often a good method to control how 
dependent are the results on the initial guesses. The length of the 
chains was chosen, after some test runs, so that it was always longer 
than the effective number of iterations required to reasonably con- 
clude that convergence has been reached. For each run, conver- 
gence was tested using graphical indicators such as the sequential 
mean and variance. A typical burn-in sequence of 5000 iterations 
was discarded at the beginning of each Markov Chain realization 
in order to ensure that the sample was indeed generated according 
to the target PPD (see Appendix [A). 

Based on test runs, we set the Gaussian components of 



the instrumental law q(0*\0 ), given by Eq. jlO) , to = 
(0.005 M G ,0.01 Gyr) and $ 2 = (0.01 M Q , 1 Gyr). At this point, 
we should note that the joint PPDs for (M, f*) show strong correla- 
tions between the two parameters, as can be seen in Fig.[T]in which 
we display them. Such a behaviour was expected from a theoretical 
point of view, but it certainly affects the efficiency of the MCMC 
algorithm, the instrumental laws used to sample the PPDs being 
symmetrical. This effect is difficult to correct manually, since the 
correlation coefficient appears to change from one combination X 
to the other and thus cannot be known a priori. 

Empirically, we notice that the composite structure of the in- 
strumental law helps us to converge if compared to simple ones, ei- 
ther Gaussian or uniform. This can indeed be seen from the accep- 
tance rates for the different components of qiO*^) from one run 
to the other. Whereas it is always low for the uniform distribution, 
it varies significantly between runs for the large and the medium 




Figure 1. Joint Posterior Probability Densities for the couple (M, t+) obtained from our MCMC simulations. Each one is labeled with the number given in the 
first column of Table|2] 
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Table 2. Inferred stellar parameters in the (M - f*) plane of the parameter space constrained using various combination, X, of the observables. First column: 
the number of the run, second column: observational constraints included in X, third column: prior used on the mass. The next columns give the estimates 
for these quantites in the form of the MP estimate and the associated credible lo" credible interval. Columns 4 and 5: estimates obtained with our MCMC 
algorithm, columns 6 and 7: estimates obtained by computing the PPD on (216 X 216) grid sampling regularly the (M - t+) plane, columns 8 to 13: same with 
(108 X 108), (32 X 32) and (16 X 16) subgrids in the (M - ?*) plane. Note that the dimension of the grid is indicative, the effective number of points used is 
lower since some models could not be computed (see text). 



MCMC Grid #1 (216 x 216) Grid #2 (108 x 108) 



Run # X 7r(M) M (M g ) 4 (Gyr) M (M G ) tl (Gyr) M (M G ) t\ (Gyr) 

1 (L,T e g) Uniform 1.121 ±0.019 5.63 ±1.14 1.120 ±0.020 5.67 ± 1.20 1.120 ±0.020 5.68 ±1.20 

2 (L,r eff ,R) Uniform 1.116 ± 0.011 5.98 ± 0.52 1.116 ± 0.012 5.99±0.53 1.116 ± 0.012 5.99±0.53 

3 (L,T eS ) Gaussian 1.107 ±0.007 6.47 ± 0.44 1.107 ±0.007 6.48 ± 0.44 1.107 ±0.007 6.48 ± 0.44 

4 (L,7" cff ,R) Gaussian 1.108±0.006 6.34 ± 0.28 1.108±0.006 6.34 ± 0.28 1.108±0.006 6.34±0.28 

5 (L, 7/ cff , <<Sv» Uniform 1.130 ±0.008 5.03 ± 0.39 1.131 ±0.008 5.00 ±0.39 1.130 ±0.008 5.03 ± 0.39 

6 (L, r cff , <Av» Uniform 1.109±0.009 6.35 ± 0.38 1.111 ±0.009 6.28 ± 0.37 1.111 ±0.009 6.28 ± 0.37 

7 (L,T eff ,R,(Sv)) Uniform 1.130 ±0.008 5.36 ±0.34 1.130 ±0.008 5.35 ± 0.35 1.130 ±0.008 5.36 ±0.35 

8 (Z.,r eff ,R,<Av» Uniform 1.103 ±0.007 6.62 ± 0.29 1.103 ±0.007 6.60 ±0.27 1.103 ±0.007 6.60 ±0.27 

9 (L,T sff ,R,(Sv)) Gaussian 1.116±0.006 5.98 ± 0.25 1.115 ±0.006 5.99 ± 0.25 1.115±O.0O6 6.00±0.25 

10 (L,r eff ,R,<Av» Gaussian 1.103 ±0.005 6.58 ± 0.20 1.104 ±0.005 6.56 ± 0.20 1.104 ±0.005 6.56 ± 0.20 



Grid #3 (32 X 32) Grid #4 (16 X 16) 

Run # X n{M) M (Ma) t\ (Gyr) M (M Q ) U (Gyr) 



1 (£,r eff ) Uniform 1.120±0.020 5.68 ± 1.19 1.120*0.020 5.70 ± 1.21 

2 (L,T eS ,R) Uniform 1.116 ± 0.012 5.99±0.53 1.1 15 ± 0.012 6.04 ± 0.55 

3 (L,T ee ) Gaussian 1.107 ±0.007 6.48 ± 0.44 1.106 ±0.007 6.45 ± 0.44 

4 (L,T eff ,R) Gaussian 1.108 ±0.006 6.35 ± 0.29 1.102 ±0.005 6.64 ± 0.23 

5 (L, r cff , <<Sv» Uniform 1.131 ±0.008 5.03 ± 0.39 1.131 ±0.008 5.02 ± 0.39 

6 (L,r eff ,<Av» Uniform 1.113 ± 0.009 6.19±0.37 1.101±3xl0~ 5 6.70 ± 1 x 10~ 3 

7 (L,r eff ,R,<5v» Uniform 1.129 ±0.008 5.38 ± 0.33 1.124 ±0.002 5.65 ±0.12 

8 (L,r cff ,R,(Av» Uniform 1.102 ±0.004 6.66 ±0.18 1.101 ±2xl0 -5 6.70±9xl0~ 4 

9 (L,r cf f,R,<(5v» Gaussian 1.115 ± 0.005 5.99±0.24 1.120±0.008 5.81 ±0.34 

10 (L,r eff ,R,<Av)) Gaussian 1.101 ±0.002 6.69 ± 0.07 1.101 ± 5 x 10~ 7 6.70 ± 2 x 10~ 5 



Gaussian. Therefore, the possibility to switch from one law to the 20 
other compensates the fact that it was not adjusted from simulation 
15 to simulation. The probabilities for law selection in Eq. l |10[ l were 
set to pi = pi = 2p3 = 0.4. Of course, this is only one possible 
solution for the choice of the instrumental law and it might not be 
the optimal one. We have merely selected it for its convenience of 
implementation and handling. 



5 5.1.2 Estimates of the parameters 

The results for our MCMC simulations are reported in Table [2] 
We considered several possible combinations for X and give them 
in the second column. The third column shows the prior we used 
on the mass. We test the effect of uniform, Eq. {4|, and Gaussian, 

10 Eq. (14) , priors on the mass. Based on the discussion from Sect. |4.4| 
we shall insist on the fact that using the latter is equivalent to make 
a strong assumption on the value of the mass compared to the non- 
informative case involving a uniform distribution. We always con- 
sidered a uniform prior on the age defined on the 2-9 Gyr domaiij^] 

15 From the MCMC simulation we directly obtain the joint PPDs 

for the parameters. This is done by simply representing the M (,) 
and t+,t= 1, . . . , T in the form of a 2-dimensional histogram (see 
Fig. [TJ. It is also possible to obtain the marginal densities for each 
parameter. These are usually computed by integrating the PPD over 



3 With the exception of run #1 in Table [2] which required a 1 Gyr lower 
limit for the uniform prior. 



the other parameters (sometimes called nuisance parameters). For 
instance, the marginal PPD for the mass is given by 

tt(M|X) = J niMj^XJdu (15) 

This integral can be approximated in the same way we did for the 
joint PPD. A quick inspection shows that they are all close to Gaus- 
sian distributions (which can be checked straightforwardly by com- 
puting the skewness and the kurtosis of the distribution). We can 
use the PM estimate and compute the Posterior Standard Devia- 
tion (PSD, square root of the Posterior Variance, noted <x) as the 
associated credible intervals. In this case, the PSD has the standard 
interpretation that a realization of the considered variable will devi- 
ate from the mean of the distribution by ±<x with probability ~68%. 
These values are given in Table [2] in column 4 for the mass and 5 
for the age. In the context of stellar physics, these statistical sum- 
maries are perhaps the most useful products of a Bayesian analysis. 
The derivation of such error bars allow us to quantify precisely how 
the various observables and their associated uncertainties impact 
the stellar parameters. In our ten runs, we explored the impact of 
various observables on the final estimates. We always included the 
effective temperature and the luminosity. We then combined them 
with other constraints. 

The first important results that stands out is the effect of in- 
cluding a Gaussian prior in Eq.[5] In this case, the estimates of the 
mass are systematically lower than when a uniform prior is used. 
Furthermore, they are always extremely close to the value M obs , 
which reflects the impact of this prior. This of course leaves open a 



10 M. Bazot, S. Bourguignon and J. Christens en-Dais gaard 




rcc/R* f* (Gyr) 

Figure 2. Left panel: Posterior probability density for the radius of the convective core of a Cen A estimated using X = (T e ff , L, (Sv)) and a uniform prior 
on the mass (run #5). The PPD is represented in the form of an histogram and as a kernel-estimated continuous density. The same plot is represented with a 
different ordinate scale in the inset. Right panel: Radius of the convective core as a function of the stellar age. The dots are the values obtained for a Cen A 
from the same MCMC simulation. The dashed lines represent the theoretical values for stars of mass, from bottom to top, M = 1.14 M Q , 1.15 M Q and 1.16 M Q . 



classical question that arises in Bayesian studies: is the used prior to 
be trusted? It is critical since including definitely has a strong 
influence on the final estimates. Given the general mass-age an- 
ticorrelation, the estimated ages are also higher when the prior is 
Gaussian. 90 

This picture is different when the prior on the mass is uni- 
form and the effects of the other observed quantities appear much 
more clearly. To understand this, we shall first briefly examine the 
effect of the various combination we used for X. The three others 
observational constraints we studied in detail are the radius, R, the 
average large separation, (Av), the average small separation, (Sv). 
They all lead to more precise estimates with respect to run #1. This 
effect is more important for the seismic observable than for the ra- 95 
dius. However, it is less than when a Gaussian prior on the mass is 
used. 

It is clear from Table [2] that two groups emerge among these 
quantities: the small separation on one hand and the radius and 
large separation on the other. The former, when included, leads to 100 
sensibly different estimates of the stellar parameters than the two 
latter. Using (dv) we find the highest mass for all our runs. Con- 
versely, the estimated age is the lowest. This seismic indicator is 
supposed to be mostly sensitive to the central regions of the star, 
which concentrate the essential of the total mass and is the locus 105 
of the nuclear burning of hydrogen, hence the best marker for stel- 
lar age. Based on these quantities, we obtain a general picture for 
a Cen A of a star with a mass — 1.13 M Q and an age roughly in 
the range ~5 - 5.3 Gyr, of course this is only true for the physics 
described in Sect. [3] and is susceptible of change for different pre- 
scriptions. 

Both the radius and the average large separation imply lower 
estimates of the mass and higher estimates for the age (with respect 
to the estimates obtained including the average small separation in 5 
X). For some cases, the estimates vary significantly, for instance 
between runs #5 and #6 or runs #7 and #8, and do not agree within 
their ltr error bars (although all runs agree if we consider 3<r error 
bars). This discrepancy is likely to be caused, as we shall see in 
Sect. |5.2| by the assumptions made on the physics of a Cen A, 10 
including in this case the setting of approximate values for X and 
a. Another source of disagreement, as discussed in Sect. [4] is that 
(Av) and R are both affected by the surface layers of the star and 
reproducing them would imply the use of a stellar evolution code 



that models accurately these regions. Finally, we also note that the 
effect of (Av) is the same as R on the PM estimates, although lesser 
in magnitude than including a Gaussian prior on the mass, whose 
observational determination might also suffer from a bias related to 
the poorly known stellar surface (Sect.|4~4|l. 



5.1.3 Structure of a Cen A 

These variations on the estimated mass and age with respect to X 
affect the general picture we can draw of the internal structure of 
a Cen A. A good example concerns the search for a convective 
core that has been a long-standing problem for this star (Turcotte & 
^hristensen-Dalsgaard|1998||Guenther & Demarque|2000| |Morel 
et al. 2000; Thevenin et al. 2002; Thoul et al. 2003; Eggenberger 
|et al.|2004||Miglio & Montalban|2005||Montalban & Miglioj 2006 1. 
It is possible to obtain as an output of an MCMC simulation the 
PPD for any function of the parameter^] We can thus estimate the 
probability of existence of a convective core P(r cc > 0) (or simi- 
larly P(m cc > 0)), with r cc the radius of the convective core (and 
m cc its mass). The distribution for the radius of the convective core 
obtained from run #5 is shown as an example in Fig. [2] We plotted 
in the form of an histogram and as a continuous density obtained 
from kernel estimation. However, we note that it is delicate to find 
an adequate representation since it is obviously discontinuous at 
r cc = 0. 

In the runs not including seismic constraints, we find a proba- 
bility of having a convective core of 10% and 2.5% for #1 and #2 
respectively. In run #3, for which the Gaussian prior on the mass 
was used, this probability becomes null. We shall note right now 
that, regardless the observation vector chosen, we never identify a 
model with a convective core when using this prior. For a uniform 
prior on the mass, we can again draw a line between estimates, de- 
pending whether (Sv) or (Av) has been included in X. In the former 
case, runs #5 and #7, we find a probability for the existence of a 
convective core of 16.5% and 1 1.2% respectively. In the latter case, 
it becomes null. As we shall see in Sect |5.2.3~| this is partly due to 



Once we have estimated the joint PPD for 0, we also have an estimate of 
the PPD of any quantity h such that h = MS). 
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Figure 3. Left column: marginal PPDs for the mass (bottom graph) and the age (top graph) of a Cen A. The observational constraints are X = (T c ff , L) and 
the prior on the mass is uniform. The histogram represent the results from the MCMC simulation. The lines represent the same marginal PPD computed from 
grids #2 (full line), #3 (dashed line) and #4 (dashed-dotted line). Right column: same as left column with X = (T e ff , L, (Sv)). 



the sensitivity of (Sv) to the innermost regions of the star, but it is 
also a bias due to our assumptions on Z and X . 

These results also allow us to derive an upper limit on the ra- 
dius and the mass of a possible convective core. If we base our- 
selves on run #5, we find r cc < 0.040 R* and m cc < 0.015 M, 
(these values are compatible with the results from run #7). Note 
also that the distribution of r cc is not Gaussian anymore. We can 
see here an obvious advantage of the Bayesian approach: if one 
wants to estimate uncertainties on such physical parameters, it is 
not possible to rely on point-estimation method such as likelihood 40 
maximum and use, for instance, Hessian inversion to derive them. 
One has to sample the entire distribution for this quantity to be in 
position to provide with robust statistics. 

45 

From Fig. [2] we can also obtain a general picture of the possi- 
ble evolutionary state of the convective core of a Cen A. It displays 
the size of the convective core as a function of time of the stel- 
lar age, the dots are the results from the MCMC simulation and 
the dashed line gives r cc = ^(f*) for M = 1.14 M e , 1.15 M Q 50 
and 1.16 M Q . It appears that the convective cores are in a growing 
phase. This mirror the fact that we tend to find younger age estimate 
when using (i5v>. The growing phase is relatively longer than the re- 
ceding one, explaining the significant probability obtained for the 
convective core existence. Similarly, the lower percentage for the 
existence of a convective core in run #7 is explained by the higher 
age and only a handful of models are being selected which are in 5 
the shorter time-scale receding phase of the convective core. 



5.1.4 Comparison with grid-based integration 

Sampling only M and t+ allowed us to compute a very dense grid 
in the region of interest of the space of parameters, Using such a 
grid, it is straightforward to integrate directly the denominator in 
the right-hand side of Eq. j5}, which simply is 

tw(X)= I 7r(0)«Sf(X|0)</0, (16) 

using classical quadrature formula (the trapezoidal rule in this 
case). Similarly, we can integrate out the nuisance parameters to 
obtain the marginal PPDs using formula of the same type as (jT3J . 

To proceed to this comparison, we computed a grid of approx- 
imate size N g = 216 x 216 in the (M - t+) plane. It was chosen 
so that the total number of computed models is close to 45000, 
which is the size of our MCMC sample when the typical bum-in 
sequence of 5000 models is removed. In fact, the total number of 
computed models is lower than N g because not all points in the 
(M - f*) plane correspond to stable solutions for our stellar evo- 
lutionary code. Stars massive enough might indeed reach the giant 
branch well before the required f*. Models producing numerical 
convergence errors during these phases of the stellar evolution were 
discarded. The effective number of computed models in our grid is 
Nf = 45533. 

Figure[3]shows a comparison between the marginal PPDs ob- 
tained with our MCMC simulations and from direct integration 
from this grid. We selected two example cases with X = (T^-, L) 
and X = (r cff ,L, (Sv)). The MCMC and grid-integrated distribu- 
tions are almost identical. This can be checked graphically or by 
comparing the first moments of the marginal densities given in Ta- 
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Table 3. Inferred stellar parameters constrained using seismic and non-seismic observables. First column: the number of the run, second column: observational 
constraints included in X, third column: prior on the mass, colums 4-8: PM estimates and associated credible intervals (the distributions are assumed Gaussian) 
for M, t-i, , Zo, a and Xq respectively. 



Run# 


X 


ji(M) 


M(M ) 


U 


(Gyr) 


Zo 


a 




12 


[L,r eff ,S,<<5v>] 


Uniform 


1.12 ±0.03 


4.1 


5 ±0.5 


0.027 ± 0.003 


1.7 ±0.1 


0.69 ± 0.02 


13 


[L,T eS ,R,(Sv)] 


Gaussian 


1.105 ±0.007 


4.f 


i ± 0.5 


0.026 ± 0.003 


1.62 ±0.08 


0.68 ± 0.02 


14 


[L,r eff ,i?,<<Sv>,<Av>] 


Gaussian 


1.102 ±0.005 


4.! 


1 ± 0.5 


0.026 ± 0.003 


1.60 ±0.10 


0.69 ± 0.02 



ble|2] They agree within a few percents at most. This was expected 
and is a good test of the self-consistency of our Bayesian approach 
that two radically different numerical strategies converge. At such a 
level of agreement, the choice of the method should be made on the 55 
basis of the computational performances. Producing dense grids of 
stellar models in 2-dimensional parameter spaces is relatively inex- 
pensive, therefore this should be the approach of choice compared 
to the MCMC. 

The fact that we are using a dense grid is critical. When con- 60 
sidering 5-dimensional parameter spaces we will not have this pos- 
sibility. Therefore, we use the computational flexibility we have in a 
lower-dimensional parameter space to illustrate the effect of sparser 
sampling on the calculation of the PPD. For this, we have selected 
subgrids in P by removing points from our original one. They have 65 
dimensions 216x216, 108x 108, 32x32 and 16x16. As mentioned 
before, some models could not be computed with ASTEC and the 
effective numbers of points are, respectively, 11 397, 936 and 245. 
We repeated the integration of the PPDs to obtain the correspond- 
ing marginal densities and their first two moments . We note that 
the original grid covers the (M -J*) plane with regular cells of di- 
mension 5.5 x 10~ 4 M Q x 37 Myr. The sparsest subgrid has cells of 70 
size 7.7 x 10~ 3 M o x 518 Myr. Even this large-sampling units offer 
very acceptable coverage of the space of parameters, in particu- 
lar in regard with commonly used grids which use steps in masP 
> 0.01 M G (see e.g. |Quirion, "Ch ristensen-Dalsgaard & Arentoft 
|20T0| >. 

The main result is that downsizing the sampling can signifi- 75 
cantly affect the outcome of the estimation process. The values of 
the posterior means and standard deviations for all selected sub- 
grids are reported in Table [2] For the first two subgrids, the re- 
sults agree fairly well with the MCMC simulations and the dense- 
grid integration. However, the results for the sparsest grid show so 
significant discrepancies with these reference estimates. The esti- 
mated posterior means can depart from the reference ones by a 
few percent (around 3-5% for the worst cases). We thus have a 
numerically-induced loss of accuracy. This down-sampling also af- 
fects the estimates of the uncertainties (which can be seen as the as 
precision on the parameters and should be distinguished from the 
accuracy). In some cases, they can be overestimated. For instance, 
in ran #1 the uncertainty on the age is 30% larger when using grid 
#4. Conversely some uncertainties can also be underestimated. This 
is significant since we see departures ranging from the MCMC and 90 
dense-grid estimates ranging from a factor ~3 (ran #5) to several 
orders of magnitudes (runs #6, #8 or #10). These very low uncer- 
tainties sometimes obtained with grid #4 simply means that only 



5 The problem of the stellar age sampling is a bit particular. Most of the 
time, these grids are database of stellar evolutionary tracks. This means 
that the stellar age is imposed by the stellar evolution code itself and not the 
user. We still needed to maintain a strict control of this parameter in order 
to compare properly with MCMC algorithm. 



few points in our sampling scheme correspond to a non-zero val- 
ues of the PPD. In these cases, the numerical integral we perform 
to compute the second moment returns unrealistic low values. This 
can be seen as a numerical error introduced by the sampling itself. 

In general, and unsurprisingly, the sharper the peaks of the 
densities, the more important this phenomenon is. This reflects 
that the sampling in the parameter space is not dense enough 
to capture properly the variations of the PPD. This in turn im- 
pacts the numerical approximation to the integrals as can be seen 
in Fig. [3] It appears clearly that the estimated densities depart 
strongly from the correctly-sampled cases (MCMC or dense grids) 
for X = (T c ff , L, (Sv)), less so for X = (T^g, L). 

Note that we chose the smallest subgrid because it has Af| ff ~ 
16 2 = 256 points. This is representative of the density of points we 
have in our 5-dimension grids. It is also noteworthy that this be- 
haviour is only indicative of the problem we might encounter when 
dealing with higher dimension. It is not possible to extrapolate fully 
the results of the two-dimensional to the five-dimensional case. 



5.2 MCMC simulations with dim(P) = 5 

5.2.1 Implementation details 

Our next step is to increase the number of parameters left free to 
vary. As noted in Sect. |5.1.T| considering only the stellar mass and 
age imply to make strong assumptions on other parameters. In par- 
ticular we needed to extrapolate the mixing-length and the initial 
hydrogen mass fraction from their solar values, which is always a 
delicate procedure that should ideally be restricted to stars lying in 
a close vicinity of the Sun in the space of parameters (such as solar 
twins see e.g., Bazot et al. 201 1). Most of the grids of stellar param- 
eters offer more than two values to sample from (for a recent exam- 
ple, see |Quirion, Christensen-Dalsgaard & Aren toft 20101. How- 
ever, increasing the dimension of P comes to a computational cost, 
namely that, all other things held constant, the number of values to 
consider in a volume element grows exponentially. 

Our goal here is two-fold. We first want to estimate M, t+, Zo, 
a and X together for a Cen A and determine the impact of relaxing 
the last three in regard with the results from Sect. |5.1| We then try to 
evaluate whether we can gain in efficiency when using the MCMC 
algorithm with respect to the grid approach in the 5-dimensional 
case. 

The physical settings of ASTEC remain the same. The ini- 
tial metallicity, the mixing-length parameter and the initial hy- 
drogen mass fraction are now left to vary. The relevant domains 
are [0.015,0.036] for Z , [1.0,3.0] for a and [0.6,0.8] for X . 
These limits were set after several trial runs. The Gaussian com- 
ponents of the instrumental law g(#*|0 w ) were set to $ t = (5 x 
10~ 3 Mq, 0.1 Gyr, 5 x 10~ 5 , 1 x 10~ 2 , 1 x 10~ 3 ) and g 2 = (2.5 x 
10~ 2 M G ,0.5 Gyr, 5 x 10~ 4 ,5 x 10~ 2 , 1 x 10~ 2 ). 

In order to limit the computational costs, we restricted the 
combinations we considered in the observable space forX. We used 
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the full set of non-seismic constraints and tested them together with 
(Sv) and (Av). For each case, we considered a Gaussian prior on 
the mass. For X = (r e fj, L, R, {Sv)) we also studied the effect of 65 
a uniform prior on the mass. The prior on the age, the metallic - 
ity, the mixing-length parameter and the initial hydrogen fraction 
are always uniform. Our Markov Chains are approximately three 
to four time longer than in the 2-dimensional case, ranging within 
150000-200000 iterations. The burn-in sequence is also longer. 70 
After testing several values for this latter, we chose to discard a 
generic 10000-iteration sequence. We also used 8 parallel chains 
for each run, increasing the computational cost but providing with 
a possibility to check the convergence of our simulations with some 
robustness. 

5.2.2 Estimates of the parameters 

First and foremost, it should be noted that the efficiency of all the 5 
strategies decreases when dim(P) increases. This has practical in- 
cidences on the subsequent statistical analysis. It becomes difficult 
to approximate correctly the PPD for five stellar parameters in a 
reasonable amount of time. This is due to two entangled factors i) 
stellar models take relatively long to compute and this can become 10 
a limiting factor when increasing the number of points to sample in 
P, ii) with increasing dimension, the shape of probability densities 
might become complex, this is expected to be even more critical 
for stellar models, which are highly non-linear in their parameters, 
hence being more difficult to sample properly from the simpler dis- 
tributions used for the instrumental law. 

A simple manifestation of these effects can first be seen in 
Fig. [5] We plotted the marginal densities of the five parameters for 
run #13. It is obvious by simply comparing these densities that one 
has to explore a much larger volume in P. This translates into a 
much higher uncertainty on the mass obtained from run #12 (com- 
pared to run #7) or on the age from run #13 (compared to run #9). 20 
Second, the departures from gaussianity become more important 
than what was observed in the 2-dimensional case. Third, one has 
to consider that these values have been obtained for chains of length 
~ 150 000. We controlled our results using multiple parallel chains 
and noticed that, in some cases, convergence does not seem to be 25 
reached. This might be due to one particular chain being stuck in 
a local maximum or not exploring efficiently enough the parame- 
ter space. Furthermore, the variance between the multiple chains 
is higher than in the 2-dimensional case (as could be seen using 
the between-chain variance B T defined in Appendix[A). Practically, 3Q 
this is reflected by the rounding of our estimates in Table[3] In more 
intuitive terms, this can be interpreted as a loss of numerical preci- 
sion between the 2 and 5-dimensional cases, run #12 being a clear 
example of this effect. 

Table|3]also offers good insight on the impact of observational 35 
data. The most striking fact is the huge effect the prior on the mass 
has. With a uniform prior, the precision on the mass decreases by 
almost an order of magnitude. We should also note the numerical 
advantage offered by the use of the Gaussian prior: by simplifying 
the target PPD it favours greatly the convergence properties of the 40 
MCMC algorithm. 

Since the mass measurements can be interpreted in terms of 
prior density, their impact is relatively straightforward to evaluate. 
The relation between the other observables and the stellar param- 
eters is much more difficult to unveil. Theoretically, it is expected 
that the average small separation should be a good constraint for 
the stellar age. To test this hypothesis without performing as many 
runs as we did in Sect. |5.1| (and also because the use of (Av) alone 



was impractical, see below), we performed another MCMC run 
(run #14) with the observable vector X = (T eB , L, R, (Sv), <Av». 
The results comparison with run #13 is particularly interesting. It 
appears clearly that the impact of the small separation is strong, and 
in fact much more constraining than the large separation. Compar- 
ing run #12 an run #13, it is clear that the mass parameter is mostly 
constrained by the prior, then it appears that the small separation 
is the observable that strongly controls the estimates of the other 
parameters. 

The use of just (Av) introduces a new problem. Studying the 
marginal PPD for t+ , it appears that most of the models lie close to 
the upper limit of the domain we defined. This means that our prior 
is very likely to be incorrect. However, when increasing the range 
for t t , the MCMC algorithm does not converge anymore. This is 
seen by the fact that our 8 chains are stuck in very different re- 
gions of the space of parameters. Note that the lack of convergence 
when sampling a larger volume of P could be identified thanks to 
our use of several parallel Markov chains. Indeed, if considered in- 
dividually, chains getting stuck in a local maximum could appear 
as the solution of our inverse problem. This shows the importance 
of a careful monitoring of the convergence of MCMC algorithm. 
This issue might be resolved with more sophisticated version of the 
MCMC algorithm {Robert & Casella|2003||Gregory|2005[[Beno^ 
[mar et al.|20"T0l|Handberg & Campante|201 \) . 

5.2.3 Structure of a Cen A 

The picture of the stellar interior we obtain from the MCMC sim- 
ulations in the 5-dimensional parameter space is sensibly different 
from those deduced from the 2-dimensional exercise. This is partic- 
ularly true when looking at the potential presence of the convective 
core. The mass prior and the atmosphere-dependent observables do 
not play such a critical role as they did in Sect |5.1.3~| in discarding 
the possibility of convection at the center of a Cen A. 

As before, we find a non-negligible probability of having 
a convective core when the only seismic observable included in 
X is the average small separation. We note that this probability 
might vary slightly from one chain to the other, but we still ob- 
tain an approximate probability of existence for the convective core 
P(m cc ,r cc + 0|X) ~ 0.40. This is much higher than the value ob- 
tain in the 2-dimensional case. Adding the average large separa- 
tion does not change this number the probability being in that case 
~42%. We can again put an upper limit on the mass and radius of a 
possible convective core. For run #12 we obtain r cc < 0.059 R* and 
m cc < 0.035 M*, for run #13 r cc < 0.055 R* and m cc < 0.030 M*. 
The difference being small, it shows that it is really the use of a pre- 
cise value of {Sv) that controls the estimated characteristics of the 
convective core. In Fig. [4] we show the posterior probability den- 
sities obtained for r cc from runs #12 and #13. We see some small 
departures from gaussianity. If we consider the distribution for non- 
null radius and mass and compute the MAP and its associated 
credible seQ we can conclude that an existing convective core in 
a Cen A would have r cc = 0.0441°;°°] R* and m cc = 0.014!° ;°°| M* 
from run #12. Again, this is very close to the values we obtain from 
run #13, r cc = 0.044+°™ 7 R* and m cc = 0.016!°°°? from. 



In this case, we could not assume that the densities were Gaussian, there- 
fore, we used the MAP as the estimate of the age and defined the (1 - rj)% 
credible set as Cj = {6k 6 P/t : tt^IX) q(t])}, with q the smallest con- 
stant such as P(C^|X) > 1 - r\, and Pj the subspace of P corresponding to 
parameter 8^. 
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Figure 4. Left panel: Posterior probability density for the radius of the convective core of a Cen A estimated using X = (T e g, L, R, (Sv) ) and a Gaussian prior 
on the mass (run #12). The PPD is represented in the form of an histogram and as a kernel-estimated continuous density. Right panel: Same as left panel for 
X = (Teg-, L, R, (6v), <Av» (run #13) 



In Fig. [5] we represent, for each individual parameters the 
range within which are observed convective cores in run #13. We 
see that they span the entire probability densities for the mass, the 
age and the mixing-length parameter. We can on the other hand 
determine lower and upper limits for, respectively, Z and X for 
the onset of convection. The resulting condition is Zq > 0.23 and 
X < 0.69. It is interesting to note that this upper limit for X does 
not correspond to the value we adopted in the 2-dimensional case. 
Whereas in Sect. |5. 1.31 we obtain some models with Xq = 0.7 and 
a convective core, it is never the case in the 5-dimensional case. 
This means that when these parameters are anticorrelated and that 
when the initial hydrogen fraction is as high as we assumed when 
we fixed it, then Z is much lower than the value of 0.027 we as- 
sumed above. Therefore, the models with convective cores that we 
observed in both case do not recover the same regions of the pa- 
rameter space P. This shows how difficult it is to make a priori 
assumptions on the stellar parameters and how misleading this can 
be in some cases. 

The general picture obtained from the 5-dimensional mod- 
elling is that the critical parameters for the existence of a convective 
core are those controlling the chemical composition. Interestingly 
enough, the precise value of the mass is not relevant for the for- 
mation of these structures, in the sense that the distribution of the 
m cc and r K for a Cen A is independent of M. This is slightly sur- 
prising since the appearance of a convective core is linked to the 
transition from the pp cycle to the CNO cycle as the main source 
of energy production at the stellar centre. On the ZAMS, when the 
temperature is high enough, the energy production rate of the latter 
(scno) overcomes the one from the former (e pp ) by several orders of 
magnitude. For a given chemical composition, £cno increases very 
rapidly with the temperature and even the modest rise of central 
temperature with mass is enough to trigger the burning of hydro- 
gen through the CNO cycle. The substantial increase in the nuclear 
energy rate that results leads in turn to an increase of V rad , hence to 
convectively unstable layers. However, in our results, we observe 
no correlation between the mass of a Cen A and the existence (or 
size or mass) of a convective core. It is noteworthy that this result 
also holds for run #12. This means that, even though the absolute 
value of the mass is important (low-mass stars will almost always 
be in the solar or M-dwarf regimes, higher-mass stars will always 
have a convective core), at the precision level obtained on this pa- 



ramete^] it is not its potential variations that will affect the devel- 
opment of the inner convective instabilities. 

We thus need to examine the influence of the other stellar pa- 
rameters on the convective core formation process. As for the mass 
we observe them for all considered values of t t and a. However, the 
situation is different for the two parameters. The mixing-length pa- 
rameter is essentially decorrelated from the mass and radius of the 
convective core. For the age, we need to look more closely to the 
distribution n(t+, m cc |X), which is represented in Fig. [6^. We can 
see that, for m cc > 0, it is skewed towards older ages. One thus can- 
not only rely on the extremum age values at which convective cores 
are observed and given in Fig. [5] since the upper limit is only rep- 
resentative of a tail of the 2-dimensional distribution tt^, m ct .|X). 
We will discuss below the implication of this preliminary remark. 

We now turn our attention towards X and Z . It must first 
be noted that these parameters are strongly anti-correlated (see 
Fig. j6jj). Therefore, in the following discussion, the effect of an 
increase in Z can be directly translated into a decrease in X . We 
should also recall that on the main sequence, since only hydrogen 
is burned and helium produced, and because we have neglected mi- 
croscopic diffusion, the initial metallicity correspond to the final 
metallicity. The main physical variables governing the energy rate 
(see e.g., Kippenhahn & Weigert 1994) are the temperature, T, the 
density p, the hydrogen mass fraction, X, and, in the case of the 
CNO cycle, the CNO-element mass fraction Xqno (when we re- 
fer to their central values, we label them with the subscript "c"). It 
is our goal here to understand how they are affected by the initial 
chemical composition of a Cen A. It can be seen in Fig. [3J; that 
ultimately Z is correlated to e c (which is £ ppc + £cno,c) and this is 
the relation we want to elucidate. 

Fig. [SJ; and |6jl clearly show that two families of models ex- 
ist for a Cen A. In the first one, the metallicity is comprised be- 
tween approximately 0.020 and 0.026. In this interval, the nuclear 
energy production rate does not increase much whereas the cen- 
tral temperature is positively correlated with the metallicity. In the 
second regime, i.e. for 0.034 > Z > 0.026, £ c increases with Z . 
On the contrary, T c does not seem to increase anymore and an ab- 



7 And even in the case of stars other than a Cen A for which the presence 
of a convective core is dubious: precisions on the mass of the order of what 
we report in run #12 are now routinely claimed from asteroseismic studies. 
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Figure 5. Left column: The histograms represent the marginal PPDs for M, t t , Zq, a and Xq obtained from run #13 (X = (T e g, L,R, (<5v)), Gaussian prior 
on the mass). The dots represent the same distribution but obtained from the integration of a grid of size ~ 16 5 . The yellow-shaded area mark the range for 
each parameter within which models with convective cores have been obtained. Right column: The white histogram are the same as in left column. The blue 
histograms are the marginal PPDs but with a reduced Markov chain. The dots represent the same marginal PPDs but obtained from the integration of a subgrid 
of ~ 8 s models 
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Figure 6. 2-dimensional marginal Posterior Probability Densities for a) the age of a Cen A, and relative mass of its convective core (note that the PPD has been 
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production rate, d) the initial metallicity and the central temperature, e) the initial metallicity and central hydrogen mass fraction, f) the central hydrogen mass 
fraction and age. 



17 



Table 4. Same as Table[3]but with estimates obtained from the grid described in Sect. |5. 2.4] 



Run# 


X 


n(M) 


A* (Mo) 


U (Gyr) 


Zo 


a 




12 


[L,T eB ,R,(6v)] 


Uniform 


1.12 ±0.03 


4.9 ±0.5 


0.026 ± 0.005 


1.6 ±0.1 


0.70 ± 0.02 


13 


[L,r cff ,R,<5v>] 


Gaussian 


1.105 ±0.007 


4.9 ± 0.5 


0.025 ± 0.005 


1.61 ±0.08 


0.70 ± 0.03 


14 


[L,r eff ,R,<av>,<Av>] 


Gaussian 


1.102 ±0.005 


5.0 ±0.5 


0.024 ± 0.004 


1.58 ±0.08 


0.70 ± 0.02 



solute maximum in the marginal PPD is reached at temperatures 
~ 1.86 x 10 7 K. 

We can interpret these two regimes as reflecting the energy 
generating cycle at play in the central layers of a Cen A. The first 
one correspond to the pp cycle. The increase of temperature with 
Z correspond to an opacity effect in the higher layer of a Cen A. 
On the one hand, enhancing the metallicity increases the opacity in 
these region^] On the other, we observe that the surface luminos- 
ity is correlated with neither Zo, e c nor T c . Therefore, to produce 
the same surface luminosity the temperature near the core needs to 
increase so that e c does too. This raise in T c has relatively modest 
impact on the energy production rate, which is expected for the pp 
cycle. 

Models of a Cen A with metallicities higher than > 0.026 
have clearly a different behaviour. We consider that they are stars 
for which the CNO cycle has already become the main source of 
nuclear energy in the central regions. The fact that the temperature 
does not increase with metallicity anymore is precisely due to the 
sensitivity of £cno to the temperature. Models with higher temper- 
atures would have had much too large luminosities for a Cen A. 
Conversely, the fact that e c becomes strongly correlated with the 
initial metallicity can be explained by an increase in ^cno.o on 
which the energy production rate depends linearly. Indeed, this 
quantity is considered to be proportional to Z and hence Z. Since 
Z CN0 can approximately double in the considered range for the 
metallicity, we are merely seeing the effect of this increase in the 
right part of the distribution in Fig.|6j;. It is extremely interesting to 
note that convective cores are present in a Cen A models almost im- 
mediately as the CNO cycle becomes dominant, the great concen- 
tration of the energy generation rate near the centre that triggering 
convection. Therefore, the CNO peaks of distributions in Fig. [5J: 



and 
Fig 



1 can be related straightforwardly to the peak for m cc > 



Finally, Fig. [6^ and]^ show that these two regimes also corre- 
spond to an age effect and we observe a clear bimodality in the dis- 
tribution for (Zq,X c ). Models with initial metallicity > 0.026 have 
a higher X c and are thus younger than those with lower metallicity. 
This is confirmed by Fig. [(J, which, besides displaying the well- 
known age-central hydrogen fraction correlation also shows that in 
the CNO mode X c decreases less rapidly with age than in the pp 
regime. This is related to the presence of a convective core in these 
models, which allows the less hydrogen-depleted regions of the star 
to act like a reservoir for the centre. We can also note that in the 
CNO regime, X c seems to be positively correlated with Z . There- 
fore, there is a slight trend for models with a higher metallicity to 
be younger. This is clearly the opposite in the pp regime. 



5.2.4 Comparison with grid-based integration 

The last step of our study was to compare the results from the 
MCMC simulations, in a similar fashion to what has been done 
in Sect. |5.1.4| to the those obtained from integration on a grid of 
models. According to our preliminary remarks, we had to compute 
a sparser grid than what we did for the 2-dimensional case. We 
chose to sample the space of parameters using 16 points in each di- 
rection. This amounts to 914 749 points in P, after discarding those 
which do not lead to a convergent solution using ASTEC. 

In Fig. [5] we represented the marginal probability densities 
as obtained from direct integration based on our grid. The results 
are fairly similar to the MCMC simulation for the mass, the age 
and the mixing-length parameter. They differ significantly for Zo 
and X . Based on our discussion from Sect. |5.1.4| we interpret 
this as a numerical effect stemming from the sparse sampling^] Ta- 
ble [4] also gives the first moment of the distribution obtained by 
the grid-integration procedure. Note that we chose to display them 
in the standard form of an estimate and the limits of an associ- 
ated symmetrical credible interval. We did with the pedagogical 
purpose in mind to show that providing such statistical summaries 
without knowledge of the underlying density might be misleading. 
Even though the posterior standard deviations are close, one can see 
from Fig.[5]that they do not recover the same reality of our state of 
knowledge. Therefore, one has to consider some of these values 
with care: the Gaussian approximation can very well be wrong, es- 
pecially when dealing with the highly non-linear stellar evolution- 
ary models. In this case of the grid-based determinations of Z and 
Xq the PM is not even the best estimator one can choose. One may 
rather want to use the MAP or the posterior median. The case for 
the metallicity is even more extreme. We see that the improper sam- 
pling leads to a multimodal distribution, and one might be tempted 
to distinguish two possible scenarios. 

These departures between the MCMC and grid-sampling ap- 
proach are, once more, of much lesser magnitude when the target 
PPD is broader. Therefore, we conclude that, since no predefine 
rule for a regular sampling of P exists, the MCMC approach is 
able, when the set of observable allows to constrain enough the 
PPD of the stellar parameters, to obtain more information on these 
densities, at an inferior computational cost. 

However, one might still be content with the results from the 
grid-based estimation and consider that they are close enough to 
what seems to be the correct PPD. Thus, the question we are now 
facing is to determine, even roughly, how much more efficient the 
MCMC could be in this 5-dimension parameter space and if this 
performance can be approached by a grid using the same number 
of model. To that effect, we selected a subgrid by rejecting one 
point out of two for each parameter. This leads us a total of N ~ 8 5 



Most likely through an increase of H~ opacity at the surface, which is 
sensitive to the metallicity. Other sources, like, bound-bound transitions or 
ionization of heavy atoms in deeper layers might also have an effect, albeit 
smaller. 



9 But of course, for real data, one cannot rule out with absolute certitude 
the possibility that the MCMC algorithm got stuck in a local maximum 
and could not sample properly the PPD. However, given the fact that our 8 
independent chains do agree, it seems highly unlikely here to be the case. 
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points sampling the volume of interest in P. In Fig. [5] (right col- 
umn), we represented the marginal PPDs obtained from the inte- 
gration of this subgrid and those from the MCMC simulation using 65 
a Markov chain of length equal to the number of points in the sub- 
grid. The result is spectacular, in the sense that the grid approach 
fails severely to approximate the densities. On the other hand, the 
reduced MCMC sample is extremely close to the large one. This 
gives a convincing illustration that, for dim(P) = 5, the stochastic 70 
sample approach offered by the MCMC methodology is much more 
efficient than the grid strategy. Our conjecture is that this remains 
valid for dim(P) > 5. 

It indeed is a well-known property of stochastic sampling 
methodologies that they tend to perform better than more sys- 75 
tematic sampling schemes for high-dimensional parameter spaces. 
However, it was our goal to verify that this was already satisfied for 
stellar parameter spaces of dimension at least 5. We also wanted to 
verify if this could translate into significant computational gains. 
This example is a very good indication that it is the case and that so 
stellar parameter estimation can benefit greatly from the use of 
MCMC strategies. 

5.3 Comparison with other studies ^ 

In this section we briefly compare our results to previous studies 
using seismic data for a Cen A (Eggenberger et al. 2004, Miglio 
|& Montalban|2005| |de Meulenaer et al.||2010| >. It is important to 
first stress the methodological differences between these previous 
studies and ours. They both rely on optimization approaches, and 90 
none of them in a Bayesian framework. This poses obvious prob- 
lems when comparing the results. First, it is not clear how fixing the 
mass compares with our use of a Gaussian prior. More importantly, 
without access to the PPDs of the parameters, it is not possible to 
assess whether or not the maximum values reported in these studies 95 
correspond to local or global extrema. We already saw in Sect. |5.2| 
that some PPDs were difficult to sample (in the 5-dimensional case, 
with X = (L, T e g, R, (Av)) and a Gaussian prior on the mass), there- 
fore such a configuration might not be uncommon. 

Different physics was also used in the stellar codes. The ma- 100 
jor point is that both studies used microscopic diffusion, which we 
have neglected. This should have an observable effect, but its nature 
remains difficult to anticipate. Gravitational settling will deplete the 
external convective zone of its heaviest elements. The initial metal - 
licity might thus be increased, or the initial hydrogen mass fraction 105 
decreased, in order to reproduce the observed Z/X. The impact of 
this additional mechanism shall be studied with care in future stud- 
ies using the MCMC approach. 

Finally, a Cen B has been modelled simultaneously in these 
studies. This introduces additional constraints on the solution to 110 
the inverse problem. We did not retain this strategy for logistical 
reasons. Modelling two stars at the same time would have implied 
significant modifications in our code. We thus preferred to focus on 
the star with the most reliable data (in particular the seismic data). 

We should also consider the effect of the value they used 115 
for the average small separation, (Sv) = 8.7 ± 0.8 yuHz used in 
these studies. To evaluate rapidly the impact this has on the es- 
timate of the age, we can use the grid described in Sect. |5.2.4| 
We obtain a relatively Gaussian marginal PPD and a lower age 
f* = 3.84 + 0.61 yuHz for X = (L,T eS ,R,(Sv)) and the Gaussian 
prior on the mass. Note that although we use the value provided 
by |de Meulenaer et al.H2010| >, (6v) = 5.8 ± 0.1 jiHz we obtained 
f* = 6.04 ±0.35 yuHz for the same X and prior. This agrees with the 
theoretical behaviour of the small separation, decreasing with age. 



We should however recall that although these grid-based are conve- 
nient because easily computed once the grid has been produced, we 
are still affected by the drawbacks we mentioned in Sect. |5.2.4| The 
marginal PPDs for the Z , a and X are obviously problematic and a 
proper MCMC simulation would be necessary (or a denser grid on 
larger scales) would be necessary to provide acceptable estimates 
of these parameters. 

Keeping in mind these remarks, we see that for the age of 
a Cen A we estimate in the 5-dimensional (runs #12 and #13) is sig- 
nificantly lower than those quoted in Miglio & Montalban ( 2005[ >. 
Their Ale model seems to be the more relevant for a comparison 
with our own results. They obtain an age of 6.7 ± 0.5 Gyr. This is 
the main departure between the two studies. The other parameters 
are all within their respective error bars (which are, as a general rule 
commensurate). However, given the complex interplay between the 
stellar parameters and other variables such as, for instance, X c , a 
change in one of them, even within the quoted error bars might 
very well lead to a very significant departure of another param- 
eter. Only an homogeneous comparison using both models, with 
the same physical prescriptions, can give us more information on 
the differences between our studies. The same kind of picture can 
be drawn from a comparison between our work and the results of 
Eggenberger et al. ( 2004"|> 

The problem of the uncertainties is also difficult to tackle. 
As a general rule, those from Eggenberger et aL| ( |2004| > seem to 
be underestimated with respect to ours. However, their method 
for estimating credible intervals is much less robust and statisti- 
cally sound than what is allowed from the study of the marginal 
PPDs of the stellar parameters. Miglio & Montalban ( 2005 ) used a 
Hessian-based approach and obtained uncertainties commensurate 
with ours. Our study also provides a confirmation of theirs since 
they need to make the implicit assumption that they are dealing with 
normally-distributed parameters before inverting their Hessian ma- 
trix. Note also that if this is not the case, this methodology is known 
to underestimate the uncertainties. 

One could also contend that the differences come from the 
physics in our models, in particular the use of the EFF equation of 
state. This is unlikely since we see in Miglio & Monta Mr7J{2005l > 
that a change of equation of state only modestly affects the final 
estimates (see for instance their models Ale and Air). Therefore, 
we suggest that most of the difference should be attributed to the 
adopted value for (6v). 

Finally, we note that gyrochronology allows to obtain some 
measurements of the age of the stars. One should note that these 
could have been used as a prior on the age. We chose not to do so, 
in particular because we preferred to constrain the age using seis- 
mic quantities rather than by the use of gyrochronology. Neverthe- 
less, we note that our estimates using the average small separations 
give an age in agreement with the results from Barnes ( 2007}. On 
the other hand, the estimates from Delorm e et al.| ( |201 1\ are much 
higher and agree more closely with our estimates using (Av). It is 
an interesting fact that the same technique provide such different 
estimates. 



6 CONCLUSION 

In this paper, we applied a Bayesian methodology to solve the in- 
verse problem that consists in estimating the stellar parameters of 
a Cen A, some observable quantities being known. To that purpose, 
we used an MCMC algorithm in its classical Metropolis-Hastings 
form. We repeated the experiment twice, first fitting the data by 
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varying only the mass and the age, and using assumptions to fix 60 
the mixing-length parameter and the initial hydrogen fraction and 
metallicity; second letting all five parameters free. 

In each case we were able to obtain the posterior probability 
densities for all considered parameters and their derived values, in 
particular the mass and radius of the convective core. We observed 65 
that different families of solutions emerged depending whether the 
observations we used are sensitive to the innermost regions of the 
star (the average small separation) or to its outer layers (all the oth- 
ers). This remarkable fact allows us in particular to reopen the case 
for a presence of a convective core in a Cen A. Expressed in proba- 70 
bilistic terms, we obtained non-negligible odds for the existence of 
such a structure, sometimes reaching over 40%. 

From a numerical standpoint, we also produced evidence that 
the MCMC approach offers a good and efficient alternative to meth- 
ods generally used in stellar physics, either optimization or grid- 75 
based Bayesian approaches. More precisely, when the dimension 
of the parameter space increases, it becomes less computationally 
expensive than the straightforward grid integration approach (as- 
suming one wishes to obtain the same statistical output). It also 
preserves the advantage of Bayesian methodologies with respect to ^ 
optimization methods in the sense that we have access to the full 
posterior probability density. Considering that it is already efficient 
and that there exist other variants of the MH algorithm that can 
improve the efficiency of our own, we deem this method as very 
interesting for the future of stellar parameter estimation. 

Finally, it seems clear that further progress on the understand- 
ing of the inner structure of a Cen A will require more precise and 
accurate data. It is in particular our hope that asteroseismology can 
help to unravel the convective core puzzle. We have already seen 
that precise values of (Sv) affect greatly our answer to this prob- 
lem. For this, one will need very long and precise radial velocity 90 
measurements, if possible continuous. The MCMC methodology 
will then be of great use to provide robust estimates of the stellar 
parameters and their uncertainties, knowing these new data. 
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APPENDIX A: MONITORING OF CONVERGENCE 
Al Algorithm monitoring 

Assessing the convergence of an MCMC algorithm to its target dis- 
tribution is a vast and complex problem. There is no clear-cut an- 
swer to this issue and most of the current methods are empirical 
ones. We present in this section the few simple indicators we used 
to monitor convergence in the present study. 

Al.l Burn-in sequence 



As seen in Algorithm |[A.1]| the initialization of the process in- 
volves the input of an initial value for the parameter. Hence, during 
a first phase, traditionally called the burn-in, the random generation 
will not be done according to the target PPD. It is then necessary to 
discard the first iterations of the Markov chain in order to yield the 
stationary regime. It is difficult to define with precision the length 
of the burn-in sequence, which obviously depends on the target dis- 
tribution and of the initial guess. 

For an MCMC sample to approximate the targeted PPD, it 
is necessary that the W are generated randoml}[^] and, hence, are 
independent of each other. 

To verify this, one can compute the autocorrelation function 
(ACF), which can be written, for a Markov Chain of length N 
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the denominator being a normalization factor. A value close to 
zero is a good indication that the values of the sequence have indeed 
been generated randomly. The behaviour of the ACF is grossly that 
of a decreasing exponential exp (—k/i p ), l p being the characteristic 
length for decorrelation (Gregory 2005 1. 

For each run, we tested burn-in sequences of size ~ l p and 
higher. Usually, discarding a number of iteration twice as large as 
l p guarantees that the sample is indeed stationarily generated. 



A1.2 Convergence testing 

As mentioned above, MCMC algorithms rely on the convergence 
properties of the Markov Chains. Therefore, the number of itera- 
tions has to be large enough in order to guarantee that the Markov 
Chain asymptotic regime is a good approximation to the target 
PPD. Assuming that an adequate burn-in sequence has been re- 
moved, the question arises of the stopping rule. 

A classical way to address the problem of convergence is to 
use empirical graphical indicators. It is particularly interesting to 
follow the evolution of the empirical cumulative average estimator 



(A2) 



where h is an arbitrary function of 6 and T the length of the Markov 
Chain. In practice, we will be mostly interested by the function 
h{0) = 6j for i = 1, . . . ,n p , so that S T approximates the Posterior 



10 Not sufficient though: one has still to determine if they are indeed gen- 
erated according to the desired PPD. 



Mean estimate |7J. At convergence, variations in S T (as a function 
of T) should be negligible. 

Depending of the dimension of the parameter space we con- 
sidered, we ran M = 3 or M = 8 chains with different initializa- 
tions. This was done in order to explore fully the space of parameter 
and eventual dependence on the initial guess. In this configuration, 
instead of S T , we consider the average over the M chains: 



(h(0)) 



(A3) 



where m indexes the number of the chain, and the associated aver- 
age variance (within-chain variance) 

.2 



(A4) 



One can also define the between-chain variance 
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(A5) 



this other indicator can be used to check convergence as explained 
in |Gelman| ( [T996) . 



